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ABSTRACT 

Perfectly matched layers are a very efficient way to absorb waves on the outer edges of media. We present a stable 
convolutional unsplit perfectly matched formulation designed for the linearized stratified Euler equations. The technique 
as applied to the Magneto-hydrodynamic (MHD) equations requires the use of a sponge, which, despite placing the 
perfectly matched status in question, is still highly efficient at absorbing outgoing waves. We study solutions of the 
equations in the backdrop of models of linearized wave propagation in the Sun. We test the numerical stability of the 
schemes by integrating the equations over a large number of wave periods. 
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1. Introduction 

The choice of appropriate boundary conditions in simu- 
lations remains a significant challenge in computational 
physics. Amidst the vast number of options that exist, a 
widely invoked construct is that of the absorbing bound- 
ary, one that is expected to relieve the computational do- 
main of outgoing waves or other structures without affect- 
ing the solution in the region of interest (for a review, see 
e.g.. IColoniusll200^ . A number of different methods exist 
to accomp l ish this, such a s an attenuating 'sponge' (e.g., 
iLuil I2OO3I : IColoniusI 2004 ). characteristics-based bound- 
ary conditions (e.g iThompson 1990}, an artificially im- 
posed supersonic outward directed advecti ve flow (e . g..lLuil 
[2003), perfectly matched layers (PMLs; iBereneeil ll994D . 



vectors perpendicular to the boundary are forced to assume 
a decaying form in the complex plane, thereby dramatically 
reducing the amplitudes of the waves in the boundary PML 
region. In the absence of d iscretization errors, the PML as 
set out by iBerenged (| 19941) is highly absorbent. After nu- 
merical discretization however, there are weak associated 
reflecti ons. An importa nt issue in the classical split formula- 
tion of[Berenger^ (I l994l ) is that the absorpt ion efficiency de- 
creases ra pidly at grazing inciden ce (e.g.. ICoUino fc MonkI 
119981 : IwTnton'fc Rappaportll2000D . 

There have been numerous advances in constructing sta- 
ble analytic continuations of wave vectors in the boundary 
region. In particular, the stable convolutional PMLs (C 
PML; e^ 



particular, the sta ble convolutional PML <s 
iRoden fc GednevI I2OOOI: iFesta fc Vilott^ I2OO5I: 



Komatitsch & Martin 
200W . 



2007 



Drossaert fc: GiannopoulosI 

contain a Butterworth filter inside the PML, 



etc. Unfortunately, a number of these methods is afflicted 
with disadvantages, a primary cause of concern being that 
of absorption efficiency. For example, the characteristics- 
based boundary conditions do not perform very well when 
the incident waves at the boundary are significantly inclined 
- and a fair degree of reflection is observed. The sponges are 
better at absorbing the waves and are relatively easy to im- 
plement; however, the reflectivity can still be substantial. 

The criterion of low reflectivity is perh aps best satisfied 
by the PML formulation, developed first by Bcrenge 3 (119941^ 
as an absorbing layer for the Maxwell equations. The cen- 
tral idea of this technique lies in performing an analytic con- 
tinuation of the wave vector into the complex plane. Wave 



thereby dramatically improving the absorption efficiency 
at grazing incidence. This formalism was adopted to 
the study of waves in anisot ropic geophysical media by 
iKomatitsch fc Martini (|2007r ) ; the numerical stability 
of sai d techniqu e was studied by Komatitsc h fc Martini 
(12001 : [M eza-Faiardo fc Papageorgiou (,2008) . Extensions 
to the poro- and visco elasti c cases have been introduce d 
by iMartin et aH (|2008[ ) and iMartin fc KomatitschI (|2009D . 
respectively. Our goal in this contribution is to develop 
a C-PML for the 3D linearized ideal MHD equations in 
stratified media. 

Astrophysical media such as stellar interiors and atmo- 
spheres may be strongly stratified and highly magnetized. 
Simulating wave propagation in such environments is of 



2 



Hanasoge et al.: Convolutional Perfectly Matched Layer 



interest because in understanding the effects of stratifica- 
tion and magnetic fields on the oscillations, we may be 
able to better interpret observations and constrain certain 
properties of the object in question. In particular, we focus 
on the Sun, a star that has been well studied and for 
which high-quality observations exist. The interior of the 
Sun is opaque to electromagnetic waves - consequently, 
only photons that arrive from very close to the surface 
(photosphere) of the Sun are visible. Helioseismology is the 
inference of the internal structure and dynamics of the Sun 
by observing the surface manifestation of i ts acoustic pul- 



satio i is (e.g., Christensen-DalsgaardI 120021 : iGizon fc BirchI 



120051: IG izon et al.ll2010ll .' Armed with accurate interaction 



theories of waves and measurements of the acoustic field at 
the surface, one can attempt to determine the sub-surface 
structure and dyna mics of various so lar features such as 
sunspots (see e.g., iGizon et al.l |2009[ ). large-scale merid- 
ional flow, interior convective length scales, etc. In order to 
construct these interaction theories, it is important to sim- 
ulate and study small amplitude (linear) wave propagation 
in a solar-like medium - filled with scatterers like sunspots 
or mean flows. The time scales of wave propagation are 
typically much smaller than the rate at which the scatterer 
itself evolves; thus one may invoke the assumption of time 
stationarity of the background medium. Numerical simula- 
tions of waves in solar- like media have been performed by 
nume r ous groups (e. g., Werne et al.l 2004 Hanasoge et al.l 



20061: IShelva g et al.l 



Cameron et al. 2007 



2006: Khomenko fc Collados 



Parchevskv fc KosovichgX 



2006; 



2007 



Hanasoge et al.l I200S ; iHanasogd 20081; ICameron et al 



2008 ) . Two groups in particular (iKhomen ko & Collados 
2006t iParchevskv fc Kosovichevll2007r currently utilize the 



classical split PML formulation in order to solve the wave 
equations in stratified media. Th ere exist drawbacks in 
these formulations. For example, IKhomenko fc ColladosI 
()2006[) have noted that there are long term instabilities 
as sociated with their method. Th e technique discussed 
in IParchevskv fc KosovichevI (|2007l ) involves the addition 
of a small arbitrary constant (see Eq. [21] and related 
discussion of their article) that could possibly be acting as 
a sponge; it is unclear if their method is perfectly matched 
after the intro duction of this constant. The m odal power 
spectrum that IParchevskv fc KosovichevI (|2007l ) display in 
Figure 9 of their article shows substantial reflected power 
from their lower boundary, atypical of perfectly matched 
formulations. 

The plan of this article is as follows: in Section 2, we 
recall the linearized ideal MHD equations and discuss the 
C-PML method as was previously developed for the seismic 
wave equations. The C-PML is then extended to the strat- 
ified Euler/MHD equations in Section 3 and results from 
numerical tests are discussed. We summarize and conclude 
in Section 4. 



2. The linearized ideal MHD equations 

As discussed in the introduction, we consider a system 
with a steady background, the wave equations are writ- 
ten as small perturbations around this state. The assump- 
tion of small amplitude wave perturbations is consistent 
with ob servations of solar and stellar oscillations (e.g., 
IChristens en-Dalsgaard 2003; Bogdan 2000). In the equa- 
tions that follow, all quantities with a subscript '0' denote 
the steady background properties while the unsubscripted 



terms represent the fluctuations around the corresponding 
property. All bold faced terms represent vectors. 



dtp 
9t(pov) 

dtp 
dtB 
V B 



-V-(pov), 
"BBo 



BqB 

47r 



S{x,y,z,t), 
-CqPoV-v - v-Vpo, 
V-(Bov-vBo), 
0. 



Bp B 

47r 



(1) 

pgo^z 

(2) 
(3) 
(4) 
(5) 



In order, these are the equations of continuity (Eq. [T]), 
momentum (Eq. [2]), adiabaticity (Eq. [5]), induction (Eq. 
U), and divergence-free field (Eq. [5]). The terms p and p 
denote density and pressure, B = {B^, By, B^) is the mag- 
netic field, go is the gravity, cq the sound speed, and v — 
{vx,Vy,Vz) is the velocity. A Cartesian coordinate system 
(x, y, z) with unit vectors (e^^. By, e^) is adopted. The iden- 
tity tensor is I = e^e^ + eyey +ezez. The momentum equa- 
tions are forced by an arbitrary spatio-temporally smooth 
source function S. All background terms po,PO)C, Bq are 
assumed to be heterogeneous, i.e., functions of (x, y, z) and 
go — goiz). Note that we have assumed the background 
medium to be non-moving, i.e., that Vg = 0. The extension 
to the Vq ^ situation requires no additional treatment at 
the boundary as long as the background velocity is zero in 
this region. 

2.1. PMLs and C-PMLs 

Consider a wave propagating in the z direction, towards 
the upper boundary. For now, ignoring the fact that the 
medium is stratified, we have Vz ~ y4e*('^^^~'^*\ where A 
is the amplitude of the wave, kz is the wavenumber along 
the vertical directio n, uj the freq uency, and i, time. The 
classical PML of Ber engeJ (Il994l) prescribes the following 
transformation for kz'. 



kz 



d ' 
1 



(6) 



where d = d{z) f {x , y , zo) > is a damping parameter 
and z = zq is the vertical coordinate corresponding to 
the entrance of the layer. We include the term f{x,y,zo) 
to account for the possibility of strong wave speed vari- 
ations in the horizontal (non-PML) directions at z — zq 
(i.e., at the entrance of the PML). Note that the damp- 
ing term applies only to fc^; i.e. the wavenumbers k^ 
and ky remain unchanged. It has been shown that, cast 
in the form of equation ([5]), the PML is unstable in a 
number of situations (for /(x, zo) = 1), in cluding in 
the p resence of mean flows (e.g., iHulfe oOl: A ppelo et al.l 
20061). and in anisotropic in edia fe.g.. iBecache et al. 20031 : 
Komatitsch fc Martini 120071 also personal communication. 



Khomenko 2009). A version of the PML with a more genera l 
transformation was introduced bv lRoden fc Gednevl (|2000D 
for Maxwell equations and adapt e d to the seismic wave 
equation by [Festa fc Vilotte ('2005VKomatitsc h fc Martini 
(200^, and Dr ossaert & Gian nopoulos (2007.) . They ap- 
plied the following relation: 



kz y kz 



(7) 
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where a — a{z) > and k = k{z) > 1 arc additional pa- 
rameters. In particular, the inclusion of a, i.e. of a filter, 
helped improve the issue of damping waves with grazing in- 
ci dence at the boundaries . We f ollow the formalism set out 
in iKomatitsch fc Martini ()2007t ) and derive the equations 
for the convolutional PML. The transformation of equa- 
tion ([7]) corresponds to a grid stretching of the relevant 
coordinate: 



(8) 



Derivatives along the vertical direction in equations ([T]) 
through ^ are now computed in terms of the stretched co- 
ordinate. In other words, the vertical derivative of a generic 
variable ip^x, y, z, t) is transformed according to: 



The term may be rewritten as: 



K {d/n + a) 



(9) 



(10) 



The first term inside the parenthesis on the right-hand 
side of equation pUj) is easily obtained by dividing the 
derivative by k. The second term is more complicated be- 
cause it involves products in the temporal Fourier space 
of two terms with frequency content; this results in a 
convolu tion in time. In o rder to compute these convo- 
lutions, 'Roden & GedneV ("20001); iFesta fc Vilotti (|2005D : 
IKomatitsch S z Martin ( 2007) use a recursion formula to up- 
date the convolution at each time step. Instead of apply- 
ing a recursive relation, here we use a differential equation 
to recover xix,y,z,uj) = sdzip (als o done in the auxiliary 
differential equatio n formulation of iGednev &: Zhao! 120101 : 
iMartin et aLll2010D . where 



s{x,y,z,uj) = - 



d/K^ 



{d/ K + a) — iuj 



(11) 



Let the time history of x such that x(a;, y, z,t < 0) = 0. 
In addition, we also require dz4' = for all t < 0. Then, 
the following differential equation for x 



dtX = ndzil^ 



- + a X, 

K 



(12) 



ensures that x h^-s the desired frequency response, i.e. 

d/K^ 



X{x,y,z,uj) 



{d/n -j- a) ~ iuj 



(13) 



The differential equations for the C-PML may now be 
written as: 



dtp 
dt (pov) 



= V 



dtp = 
dtB 



V;i-(pov) - PodiVz 
An 47r 
An An 

pgo^z - o-pov + s, 

CoPoVh-v - c^podiV; 
Vh ■ (Bov - vBo) -I- 



Vzdipo, 
Bo B 

P 



(14) 



An 
Bp B 

An 

(15) 

- vVhPo - v^dspo, (16) 
dz {Bozv - VzBo) , (17) 



where V/i — e^dx + ^ydy represents the derivatives in the 
non C-PML directions, and ga^ dzPo, d^po are the modified 
background gravity, pressure and density gradients in the 
absorption region respectively. A sponge-like decay term 
— crpov, where a — a{x,y,z) is introduced in the momen- 
tum equation in order to stabilize the system. It is impor- 
tant to note that in the absorbing layer, the solution is 
non-physical and therefore discarded. Thus any modifica- 
tions applied to the equations, in relation to issues such as 
altering the stratification or not enforcing V • B = 0, is 
justifiable as long as the solution in the interior domain is 
unaffected and the formulation is stable. The equations (|T^ 
through (fT7|) in conjunction with equation (|T^ provide the 
following prescription: 



dtp = - V;j-(pov) - po 



dt^ = --^dzVz - (-+a]^, 



diPo 



/2 

a/ K. 
d/n + a 



Vzdipo, (18) 
(19) 



dzPo, 



(20) 



9t(pov) = V/i 
+ -dz 



[BBo ^ 


BoB 


An 


An 


BzBq 


1 BozB 


An 


An 



Bp B 

An 
, Bo -B 



i - Pgo^z - crpoV + S, 



dti = 



go 



;d. 
d 

K 



BzBq BqzB 



An 



An 



An 

Bp B 

An 



d/ K + a 

dtp = -CgpoV/i-v - c^po 

- v-VhPo - VzdiPQ, 

dsPa = 77 — ■ — dzPo, 
d/ K + a 

dtB = V,, • (Bov - vBo) 

+ -dz {BozV - VzBq) + r 

K 

dtri = — '^dz (Boz^ ~ VzBo) - 



-dzVz 
K 



a q. 



(21) 



(22) 
(23) 

(24) 
(25) 

(26) 
(27) 



Note that ^, ry are vector memory variables required for 
the calculation of the convolutional part of equation (flUI) 
and that ij — (ryj., ryj,, 0). These are additional variables 
are needed in order to solve the auxiliary differential equa- 
tions (fro|) . (1 ^ . and (P7|) . 

Why is the stratification altered in the boundary re- 
gion? Consider for example, the derivative term d^po. The 
differential equation (fl^ when applied to the steady dzPo 
term leads to: 



diPo = -dzPo + X 

K 

dtX = — '^dzPo - ( - + a 1 X- 



(28) 
(29) 
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Since we are dealing with a steady background medium, we 
drop the time derivative in equation (|29p and obtain 

dzPn = 9zPo- (30) 

a/ K + a 

Note that we can pursue a similar formalism for the 
gradient term 9jpo- To maintain hydrostatic stability in 
the absorption layer, we alter gravity by prefixing it by 
the same term, -jj-j^^- And hence the relations in equa- 
tions (gni), (El), and (P5|) . 

In the presence of magnetic field, three wave branches 
appear: magneto-sonic slow, fast, and Alfven (the distinc- 
tions between these modes apply only if the magnetic 
and hydrodynamic pressures are not in equipartition) . The 
physics of the propagation of these waves is intricate and 
inher ently highly anisotropic (e.g., iGoedbloed fc PoedtsI 
I2004D . Alfven waves are transverse incompressible oscil- 
lations whereas the slow and fast modes are compressive 
magnetically guided/influenced oscillations. The slow mode 
dispersion relation exhibits certain pathologies wherein for 
certain wavenumbers, the phase speed vanishes while the 
group speed becomes very small (e.g., see pp. 195-214 of 
IGoedbloed fc PoedtsI |2004 ). This strong anisotropy desta- 
bilizes the calculation in the entry regions of the boundary 
layer. Thus to offset the possibility of this instability, we 
add the sponge term — crpov to the momentum equations 
in regions of the magnetic field. 

The choice of the parameters k, d, a is fairly central 
to creatin g an efficient absorpti o n bou ndary formulation. 
Following iKomatitsch fc Martini ()2007D . we set a = tt/q, 
where /o is a characteristic frequency of the waves. The a 
decays linearly to zero over the absorption region, being 
maximum at the start of the layer and falling to zero at the 
boundary. The damping function d = f{x, y, zq) {z/L)^ , is 
zero at the boundary layer entry point and attains its max- 
imum value at the boundary. We find that for the Euler 
equations, the following choice for / results in efficient ab- 
sorption (e.g., Collino & Tsogka 2001): 

N +l_ 

f{x, y, zo) = ^^^l°gio ^c, (31) 

where L is the length of the layer, Rc is the amount of 
reflection tolerated, and c is the characteristic sound speed 
in the C-PML layer. Note that c is not a function of z. For 
the ideal MHD equations, we use: 

N +1 

f{x,y,zo)^ — — logio (32) 



where, c^,(a;, y, zq) = a/ co{x, y, zq)^ + ca{x, y, zq)^ is the 
characteristic propagation speed of the fastest waves in the 
absorption layer (set for a specific vertical location z = zq) 
and Rc,N,L retain the same meaning as for the Euler 
equations. Wave propagation speeds can vary substan- 
tially from strongly magnetized regions to non-magnetic 
regions; for this reason, we allow Cw to vary horizontally, 
i.e., Cw = Cw{x, y, zq). The sound speed is denoted by c and 
the Alfven speed by ca = l|B||/-v/47rpn. The k te r m acts to 
damp evanescent waves (|Roden fc GednevI 120001 : iBerengeij 
1200211 : although there are evanescent waves in our simula- 
tions, we find that instabilities are set into motion when 
K is allowed to vary from 1 at the absorption layer en- 
try to 8 at the boundary. For values of k less than this. 



the improvement in absorption efficiency is not very no- 
ticeable. We therefore set k = 1. The damping sponge is 
cr = -[(7V-hl)cA/L](z/L)^logio Rc. The perfectly matched 
status of the MHD formulation comes into question with 
the introduction of the sponge term. Note that because 
a oc CA the technique is perfectly matched for the strat- 
ified Euler equations, i.e. for the ca = case. 

3. Numerical results 

3.1. Methods 

The simulation code employed in these calcul ations has 
been discussed in var ious previous pu blications (jHanasogd 
120071: iHanasQge et al.|[2008t lHanasogell2008D . The vahdation 
and verification of the numerical aspects of the code may 
also be found in these articles. In summary, to capture vari- 
ations in the vertical direction, we apply sixth order, low 
dissip ation and dispersion compact finite differences (jLeld 
|1992[ ). A non-uniformly spaced grid is used because the 
eigenfunctions change much more gradually deeper within 
the interior; the vertical grid spacing varies from several 
hundred kilometers deeper within to tens of kilometers 
in the near-surface layers. Because the dissipation of the 
scheme is low, spectral blocking occurs along the vertical 
directions (energy buildup at the highest wavenumbers). 
This is a direct consequence of the strong vertical strati- 
fication (terms like c,po,po vary strongly with z) and the 
non-uniform grid spacing; in order to avoid instabilities, 
we regularly de-alias the variables along th e z direction ac- 
cordin g to the technique described in Hanasoge fc Duvall 
()2007D . Zero Dirichlet upper and lower boundary condi- 
tions are enforced at the outer edge of the C-PMLs for all 
the variables. A Fourier spectral method is applied in the 
horizontal directions in order to compute derivatives. Time 
stepping is achieved through the app lication of a se cond- 
order optimized Runge-Kutta scheme (jHu et al.lll996t ). The 
horizontal boundaries are chosen to be periodic while the 
vertical boundaries are required to be absorbent. We study 
this case in detail below; there is no loss of generality how- 
ever, since the method can be extended in order to render 
the horizontal boundaries absorbent as well. 

An inspection of equations (ITS)) to (|27p reveals that 6 
memory variables are needed: for Vz, the three (vector) 
Lorentz force components, and two in the induction equa- 
tion. The memory and computational overheads are rela- 
tively small: (1) the boundary layer works well with 8-10 
grid points; thus the six memory variables need only be 
stored over a small number of points, and (2) the additional 
computations are limited to a small number of multiplica- 
tions and additions, as required for the evolution of the 
convolution equations p^ . (j^ . and (P7)) . 

3.2. Waves in a non-magnetic stratified fluid 

We first demonstrate the ability of the C-PMLs at absorb- 
ing outgoing waves. The background medium is a stratified 
polytrope with index m = 2.15; the vertical extent of the 
computational domain is such that approximately 2.6 scale 
heights in density and 3.72 scale heights in pressure fill the 
domain (i.e. the pressure at the bottom is e^'^^ = 41.5 
times the value at the top). The stratification properties 
of the polytrope are described in appendix lA.ll Waves 
are locally excited in a small region located at a distance 
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Fig. 1. Snapshots of the normahzed wave velocity ^/PqVz at 
four time instants, non-dimensionahzed to span the range 
±10. The grey scale in all the panels is identical. The ve- 
locities at < = 18 minutes are barely discernible. 



of 18 Mm above the bottom boundary (the full vertical 
length is 68 Mm). For the C-PML parameters, we choose 
iV = 2, i?c = 0.1%, /o = 0.005 Hz in equations ^ and ^ 
for these calculations. Both the upper and lower C-PMLs 
are 10 grid points thick. In Figure [TJ we show snapshots of 
the wavefronts at four instants of time. The scaling in all 
the plots is identical. It is seen that the waves are almost 
completely absorbed by the upper boundary. In order to 
study the absorption more quantitatively, we plot the fol- 
lowing energy invariant (summed over the entire grid) as a 
function of time in Figure H fe.g.. lBogdan et al.lll996D : 



^Po||v| 



(33) 



In order to test the stability of the C-PMLs, we integrate 
the hydrodynamic wave equations for around 300 wave pe- 
riods (peak frequency ^ 6 mHz; 12 hour time integration). 
In this case, the setup is more realistic; the computational 
domain now straddles the solar photosphere, a region where 
the density and pressure drop exponentially with height. 
The waves propagate in a domain that contains 21 scale 
heights in density, representing a contrast of 1.3 billion be- 
tween the bottom and top. For details of the stratification, 
see appendix lA. 21 

We solve these equations in a computational domain 
of size 200 x 200 x 35 Mm'^, (horizontal sides followed 
by vertical length) (1 Mm — lO^m), discretized using 
256 X 256 X 300 grid points. Vertically, the box extends 
from approximately 34 Mm below the solar photosphere 
to 1 Mm above. Standing waves or normal modes are natu- 
rally created by the stratification of the medium; the modes 
are trapped between the surface and a specific depth be- 
low that depends on the frequency and wavenumber of the 
mode. A consequence of this trapping is that the mode has 
an evanescent tail that extends from the surface into the 
overlying atmosphere. We therefore place the upper com- 
putational boundary high enough above the surface so that 
the normal mode evanescent tails are largely unaffected by 
it. 

An additional aspect is that of the wave excitation 
mechanism: in the Sun, waves are created by the ac- 
tion of the strong ne ar-photospheric turbulence (e.g., 
[St em fc Nordlundl[2000l) . To mimic this process, we add 
a horizontally homogeneous vertically localized determin- 




20 40 60 80 100 
Time (iriiiiuteH) 

Fig. 2. Time evolution of the normalized modal energies as 
computed with equation ([55)1 . The initial drop in energy 
is due to the arrival of waves at the lower C-PML from 
the source. The energy of these downward propagating 
waves is plotted as the dot-dash line. The second drop in the 
energy corresponds to the first arrivals at the upper C-PML 
(displayed with dashed lines). Modes with large horizontal 
wave numbers and weakly reflected waves arrive gradually 
at the boundaries at later times, resulting in an extended 
energy decay. Both upward and downward propagating 
wave packets are efficiently absorbed by the CPMLs. 



istic forcing function (S of the momentum equations). 
This source term is only vertically directed, i.e., S = 
S{x,y,t)5{z — Ze)ez, where Ze = —50 km. In order to gen- 
erate S, we use a Gaussian random number generator to 
populate its Fourier domain and the resultant function is 
multiplied by the average frequency power spectrum of the 
Sun, modeled as a Gaussian with peak frequency at 3 mHz 
and a full width of 1 mH z. The inverse transform of this 
funct ion is S{x,y,t) (e.s.. iHanasoge et al.ll2006t iHanasogd 
l2008h . 

The calculation is stable over the entire time integra- 
tion. The wavenumber-frequency power spectrum of the 
vertical velocity component Vz at a height oi z = 200 km 
above the photosphere is plotted in Figured The horizontal 
axis represents the non-dimensional horizontal wavenumber 
kfiRQ, where Rq is the solar radius and kh is the horizontal 
wavenumber, while the vertical axis is frequency in mHz. 

In order to study the impact of these absorbing lay- 
ers, we perform simulations with C-PMLs placed adjacent 
to the upper and lower boundaries (Figure ^ and damp- 
ing sponges (Figure |4]) . A demonstration of the efficacy of 
the C-PML is the absence of artifacts related to the lower 
boundary in the modal power spectrum in Figure [31 Weak 
reflections from the lower boundary cause the dispersion 
relation to change, manifested in the power spectrum as a 
flattening of the curvature of the ridges, as seen in FigureHl 



6 



Hanasoge et al.: Convolutional Perfectly Matched Layer 



Power spectrum 




200 400 600 800 1000 1200 1400 
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Fig. 3. Modal power spectrum of the vertical velocity com- 
ponent extracted at a height of z = 200 km above the 
photosphere, from a 12 hour long simulation with C-PMLs 
placed at the upper and lower boundaries (plotted on a lin- 
ear scale). The horizontal axis is the normalized wavenum- 
ber, the vertical is frequency, and regions of high power 
are the normal modes that appear in the calculation. The 
symbols overplotted on the power contours are the theoreti- 
cally expected values of the resonant frequencies, computed 
using the MATLAB boundary value problem solver bvp4c. 



Power spectrum 
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Fig. 4. Modal power spectrum extracted of the vertical ve- 
locity component extracted at a height oi z = 200 km above 
the photosphere, from a 12 hour long simulation with the 
boundaries lined by sponges (plotted on a linear scale). 
The weak reflections from the lower boundary cause the 
dispersion relation to change in curvature at all points for 
which i//(fc?ii?©) exceeds a certain threshold. The symbols 
are the same as those in Figure |3l note that there are differ- 
ences between the C-PML and sponge layered simulations 
at high frequencies. This is because these high-frequency 
waves are sensitive to the upper boundary. 



3.3. Stratified MHD fluid 

We first perform a test akin to that of Figure [5] A con- 
stant inclined magnetic field is embedded in the polytrope 
of appendix lA.ll the strength of the field is such that the 
Alfven speed is approximately four times the sound speed 




40 60 

Time 



Fig. 5. Time evolution of the normalized MHD wave energy 
as computed using equations summed over the entire 
grid. Because of the strong dispersive nature of the waves 
and large variations in the propagation speeds between the 
fast, slow, and Alfven waves, some modes take a long time 
to reach the boundaries, resulting in a relatively extended 
energy decay process compared to Figure^ The initial drop 
at t ~ 20 — 40 min corresponds to the first arrivals of waves 
at the upper and lower absorption layers. The secondary 
decline may be associated with the arrivals of the down- 
ward propagating dispersing slow and Alfven modes at the 
absorption layer adjacent to the lower boundary. 



at the upper boundary. Waves are excited in the vicinity 
of the upper boundary by horizontally shaking the field 
lines. The fast and Alfven waves reach the upper boundary 
first but the slow modes follow behind, which makes the 
demonstration of the absorptive properties of the bound- 
ary formulation less evident than in Figure[2l Furthermore, 
constructing energy invariants for MHD waves in stratified 
media is not a trivial affair due to the fairly complicated 
Lorentz force terms. A number of autho rs have studi ed this 
probl em (e.g.. lBrav fc Lou gh headlll974tlParker.i,1979jLerovl 
Il985[ ) and have arrived at differing versions of an invariant 
(personal communication, P. S. Cally 2009); we use the fol- 
lowing appro ximate form (kinetic -|- t hermal -|- magnetic 
energies; e.g., iBrav fc Loughheadlll974h : 



IBI 



2P0IIVI 



27P0 



Stt 



(34) 



We plot the time history of the energy summed over the 
entire computational domain in Figure [SJ 

Our final test consists of studying the long-term stabil- 
ity of the method numerically. A magnetic flux tube (e.g., 
iMoradi et al.ll2009t ) is embedded in the stratified polytrope 
(appendix IA.2|) : waves are excited in the non- magnetic 
regions and allowed to propagate through the flux tube, 
whose geometry is shown in Figure [S] The system is sta- 
ble up to 12 hours of time integration. The modal power 
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Fig. 6. The 'sunspot' like magnetic field configuration used 
in the long-term integration numerical test. On the left 
panel we plot the ratio of the local Alfven to the sound 
speed as a function of x and z. The magnetic field becomes 
dynamically unimportant in the deeper layers because the 
background hydrostatic pressure increases with relative ra- 
pidity. The right panel displays a horizontal cut (at the solar 
photosphere) through the fast mode speed {\/ c\ + Cq ) . 



spectrum from this data looks almost identical to that in 
Figure [3] (with the exception of small frequency shifts) and 
is hence not shown. The maintenance of V • B = numer - 
ically is an essential task of the scheme (e.g., iTothI [2000[) . 
Discretization errors are one of the primary causes for non- 
zero V • B = 0; additionally, in this case, the boundary 
formulation in non V • B = conserving. Consequently, 
it is important to test if the error grows with time and to 
quantify its level. In order to measure the normalized er- 
rors in maintaining V • B = 0, we compute the following 
quantity: 



£2 



i::dz\\v -Bi 



(35) 



where zi , Z2 are the limits of the physically-relevant domain 
(i.e., the non- absorbing regions) and the following definition 
for the L2 norm is applied: 



L2[f{x)]= 



(36) 



where the sum is over all the grid points Xi in the computa- 
tional domain. We display the time evolution of rig on the 
left panel of Figure [71 One can see that rie is nearly con- 
stant, perhaps showing a slight decrease with time. On the 
right panel we show the time history of the numerator of 
Tie (Eq. |35|). The gradual increase in the magnetic energy 
is due to constant input wave excitation via the S term in 
the vertical momentum equation. The growth rates of the 
magnetic energy demonstrate the numerical stability of the 
system over the period of time integration. 

4. Conclusions and Future Work 

We have developed a Convolutional Perfectly Matched 
Layer (C-PML) for the stratified linearized Euler equations 
and a highly efficient absorption layer for the ideal MHD 
equations. This boundary formulation is quite useful for the 
calculation of wave propagation in astrophysical plasmas, 
where stratification and magnetic fields abound. The ab- 
sorption layer for MHD waves is a slightly altered version 
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Fig. 7. Plotted on the left panel is the error in maintaining 
||V • B|| = as measured by rie, defined in equation (j35|) . 
On the right panel is plotted the L2 norm of f^^ dz 1 1 V • B| | 
(numerator of Eq. [35]) as a function of time. 



of th e convolutional formulation of fe.g.. lRoden fc GednevI 
l2000f ). requiring an extra sponge- like term in order to stabi- 
lize the system. While the boundary methods as discussed 
here are perfectly matched for the stratified Euler equa- 
tions, the inclusion of the sponge term casts some doubt 
as to whether the same could be said of the MHD equa- 
tions. Simulations with the absorption layers for the strat- 
ified MHD and Euler equations were performed and found 
to be stable over a time length of 300 wave periods. With 
increas ingly sophistica ted boundary formulations, such as 
that of' Hu et al.l (|2008[) . who developed the PML for turbu- 
lent flows, we may in the near future be able to extend these 
techniques to fully non-linear MHD turbulence/dynamo 
calculations. 
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Appendix A: Stratification 

A.l. Only polytrope 

For the absorption tests in Sections 13 . 21 and [3.31 we use the 
following polytropic stratification prescription: 



Po{z) 
Po{z) 
90 



Ppoly 



1 



Z 



Ppoly 



1-A 



m+1 



1 Ppoly 



Zf Ppoly 
1 



(A.1) 



Fi = 1 + -, 



co{z) 



Ppoly 



z 



(A.2) 



We set Ppoiy = 1.178 x 10^ dynes cm^^,/9poiy = 3.093 x 
10"'^ g cm-^, Zf = -0.450 Mm, and m = 2.150. Note also 
that z is the height, i. e. the atmospheric layers are de- 
scribed by z > and vice versa; 2; = is the fiducial sur- 
face. The polytrope is truncated at z = —31.77 Mm (upper 
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boundary) and z = —104.40 Mm (lower boundary) respec- 
tively. 

A.2. Polytrope + isothermal layer 

Waves propagating toward the surface in the Sun are re- 
flected by a sharp fluctuation in the b ackground density 
gradient (e.s.. IChristensen-Dalsgaardl2003i) . This reflection 
zone is located at a height of Zr ~ —0.050 Mm below the 
surface. We mimic this by patching the polytropic strati- 
fication of section lA.ll with an overlying isothermal layer. 
This patching results in a fluctuation in the density gra- 
dient, leading to an acoustic cut-off frequency of approxi- 
mately 5.4 mHz, similar to the solar value. For z < zm, we 
use the prescription of section lA.il For z > Zr, we use the 
following equations: 



Po{z) 



H 



Piso exp 




H 


Piso exp 




H 


Ppoly ^1 


Zr\ 
^f) 


Ppoly (^1 


Zr\ 


Piso 




50Piso 





(A.3) 



The relations for p-^so and piso arise from the requirement of 
continuity of the pressure and density at the matching point 
between the polytropic and isothermal layers. The relation 
for is a consequence of enforcing hydrostatic balance. 
This model is truncated at z = —34 Mm (lower boundary) 
and z = +1 Mm (upper boundary). 
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